clear all;

lambda_SOC=0.1*0.5;

spin_Delta=0.0;
spin_nvec=[0,0,1];
spin_nvec=spin_nvec/norm(spin_nvec);

sigx=[0,1;1,0];
sigy=[0,-i;i,0];
sigz=diag([1,-1]);
sig0=eye(2);

gen_SOC_term=@(bdvec) sqrt(2)*(i)*lambda_SOC*(bdvec(1)*sigx+bdvec(2)*sigy+bdvec(3)*sigz);


a1=[2.7508   -1.5882    4.4921];
a2=[0    3.1764    4.4921];
a3=[-2.7508   -1.5882    4.4921];

Rh_pos=[0.125000 0.125000 0.625000 ;
    0.625000 0.125000 0.125000 ;
    0.125000 0.625000 0.125000 ;
    0.125000 0.125000 0.125000];

Rh_pos=Rh_pos-0.125;

Rh_posC=Rh_pos;
Rh_posC(:,1)=Rh_pos(:,1)*a1(1)+Rh_pos(:,2)*a2(1)+Rh_pos(:,3)*a3(1);
Rh_posC(:,2)=Rh_pos(:,1)*a1(2)+Rh_pos(:,2)*a2(2)+Rh_pos(:,3)*a3(2);
Rh_posC(:,3)=Rh_pos(:,1)*a1(3)+Rh_pos(:,2)*a2(3)+Rh_pos(:,3)*a3(3);



Rh_posC





%%

tetra1_posC=Rh_posC
tetra1_C=mean(tetra1_posC,1)

tetra2_posC=-Rh_posC
tetra2_C=mean(tetra2_posC,1)


newx_axis=zeros(4,3);
newy_axis=zeros(4,3);
newz_axis=zeros(4,3);
for ind=1:4
    zvec=tetra1_posC(ind,:)-tetra1_C;
    newz_axis(ind,:)=-zvec/norm(zvec);
end
newz_axis

newx_axis(4,:)=[1,0,0];
newx_axis(3,:)=[1,0,0];
newx_axis(2,:)=[-0.5,-0.5*sqrt(3),0];
newx_axis(1,:)=[-0.5,0.5*sqrt(3),0];


for ind=1:4
    xvec=newx_axis(ind,:);
    zvec=newz_axis(ind,:);
    newy_axis(ind,:)=cross(xvec,zvec);
end


dot(newy_axis(2,:),newz_axis(2,:))


%%
Rh_pos(2,:)+(Rh_pos(2,:)-Rh_pos(1,:))+[-1,0,1] %
Rh_pos(2,:)+(Rh_pos(2,:)-Rh_pos(3,:))+[-1,1,0]
Rh_pos(2,:)+(Rh_pos(2,:)-Rh_pos(4,:))+[-1,0,0]

%%
Rh_pos(1,:)+(Rh_pos(1,:)-Rh_pos(2,:))+[1,0,-1]
Rh_pos(1,:)+(Rh_pos(1,:)-Rh_pos(3,:))+[0,1,-1]
Rh_pos(1,:)+(Rh_pos(1,:)-Rh_pos(4,:))+[0,0,-1]
%%
Rh_pos(3,:)+(Rh_pos(3,:)-Rh_pos(1,:))+[0,-1,1]
Rh_pos(3,:)+(Rh_pos(3,:)-Rh_pos(2,:))+[1,-1,0]
Rh_pos(3,:)+(Rh_pos(3,:)-Rh_pos(4,:))+[0,-1,0]
%%
Rh_pos(4,:)+(Rh_pos(4,:)-Rh_pos(1,:))+[0,0,1]
Rh_pos(4,:)+(Rh_pos(4,:)-Rh_pos(2,:))+[1,0,0]
Rh_pos(4,:)+(Rh_pos(4,:)-Rh_pos(3,:))+[0,1,0]
%%

tetra1=Rh_posC;
tetra1C=mean(tetra1,1);

tetra2=Rh_posC;
tetra2(1,:)=tetra2(1,:)-a3;
tetra2(2,:)=tetra2(2,:)-a1;
tetra2(3,:)=tetra2(3,:)-a2;

tetra1
tetra1C=mean(tetra1,1)
tetra2
tetra2C=mean(tetra2,1)


%%
% 
% figure(1);
% clf;
% 
% 
% hold on;
% 
% linecc= [.5 .5 .5];
% %plot_tetra(Rh_posC,0*a1,linecc)
% % plot_tetra_v4(Rh_posC,a1,linecc,2)
% % plot_tetra_v4(Rh_posC,a2,linecc,2)
% % plot_tetra_v4(Rh_posC,a3,linecc,2)
% % 
% % %plot_tetra(-Rh_posC,0*a1,'r')
% % %plot_tetra_v2(-Rh_posC,a1,linecc,2)
% % %plot_tetra_v2(-Rh_posC,a2,linecc,2)
% % %plot_tetra_v2(-Rh_posC,a3,linecc,2)
% % 
% % plot_tetra_v4(-Rh_posC,a1+a2,linecc,2)
% % plot_tetra_v4(-Rh_posC,a2+a3,linecc,2)
% % plot_tetra_v4(-Rh_posC,a1+a3,linecc,2)
% 
% 
% 
% hexpts1=[];
% hexpts1(1,:)=Rh_posC(3,:)+a1;
% hexpts1(2,:)=Rh_posC(2,:)+a2;
% hexpts1(3,:)=Rh_posC(1,:)+a2;
% 
% hexpts1(4,:)=Rh_posC(3,:)+a3;
% hexpts1(5,:)=Rh_posC(2,:)+a3;
% hexpts1(6,:)=Rh_posC(1,:)+a1;
% 
% 
% hexpts2=[];
% hexpts2(1,:)=Rh_posC(4,:)+a1;
% hexpts2(2,:)=Rh_posC(1,:)+a1;
% hexpts2(3,:)=Rh_posC(2,:)+a3;
% hexpts2(4,:)=Rh_posC(4,:)+a3;
% hexpts2(5,:)=Rh_posC(1,:);
% hexpts2(6,:)=Rh_posC(2,:);
% 
% 
% 
% hexpts3=[];
% hexpts3(1,:)=Rh_posC(4,:)+a2;
% hexpts3(2,:)=Rh_posC(2,:)+a2;
% hexpts3(3,:)=Rh_posC(3,:)+a1;
% hexpts3(4,:)=Rh_posC(4,:)+a1;
% hexpts3(5,:)=Rh_posC(2,:);
% hexpts3(6,:)=Rh_posC(3,:);
% 
% 
% 
% hexpts4=[];
% hexpts4(1,:)=Rh_posC(4,:)+a3;
% hexpts4(2,:)=Rh_posC(3,:)+a3;
% hexpts4(3,:)=Rh_posC(1,:)+a2;
% hexpts4(4,:)=Rh_posC(4,:)+a2;
% hexpts4(5,:)=Rh_posC(3,:);
% hexpts4(6,:)=Rh_posC(1,:);
% 
% hexpts1(7,:)=hexpts1(1,:);
% for indh=1:6
%     pt1=hexpts1(indh,:);
%     pt2=hexpts1(indh+1,:);
%     line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',[0,0,1],'LineWidth',7)
% 
% end
% %pp=fill3(hexpts1(:,1),hexpts1(:,2),hexpts1(:,3),[185,204,255]/255);pp.FaceAlpha=0.7;
% pp=fill3(hexpts1(:,1),hexpts1(:,2),hexpts1(:,3),[0.4,0.4,1]);pp.FaceAlpha=0.5;
% 
% 
% coorcenter=[-4.,-3.,hexpts1(1,3)-1];
% coordinatelen=1.5;
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),coordinatelen,0,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,coordinatelen,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,0,coordinatelen,'MaxHeadSize',15,'Color','k','LineWidth',7)
% 
% 
% mirror_n1=(hexpts1(6,:)+hexpts1(1,:))/2-(hexpts1(3,:)+hexpts1(4,:))/2;
% mirror_n1=mirror_n1/norm(mirror_n1);
% mirror_n2=[0,0,1];
% 
% pltm1=4.;
% pltm2=1.0;
% pltmirror=zeros(4,3);
% pltmirror(1,:)=-pltm1*mirror_n1-pltm2*mirror_n2;
% pltmirror(2,:)=pltm1*mirror_n1-pltm2*mirror_n2;
% pltmirror(3,:)=pltm1*mirror_n1+pltm2*mirror_n2;
% pltmirror(4,:)=-pltm1*mirror_n1+pltm2*mirror_n2;
% pltmirror(:,3)=pltmirror(:,3)+hexpts1(1,3);
% %pp=fill3(pltmirror(:,1),pltmirror(:,2),pltmirror(:,3),'g');pp.FaceAlpha=0.3;
% 
% pt1=-pltm1*mirror_n1+[0,0,hexpts1(1,3)];
% pt2=pltm1*mirror_n1+[0,0,hexpts1(1,3)];
% pp=line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',[0,1,0,0.3],'LineWidth',7)
% %pp.FaceAlpha=0.3;
% 
% mirror_n1=(hexpts1(2,:)+hexpts1(1,:))/2-(hexpts1(5,:)+hexpts1(4,:))/2;
% mirror_n1=mirror_n1/norm(mirror_n1);
% mirror_n2=[0,0,1];
% 
% pltm1=4.;
% pltm2=1.0;
% pltmirror=zeros(4,3);
% pltmirror(1,:)=-pltm1*mirror_n1-pltm2*mirror_n2;
% pltmirror(2,:)=pltm1*mirror_n1-pltm2*mirror_n2;
% pltmirror(3,:)=pltm1*mirror_n1+pltm2*mirror_n2;
% pltmirror(4,:)=-pltm1*mirror_n1+pltm2*mirror_n2;
% pltmirror(:,3)=pltmirror(:,3)+hexpts1(1,3);
% %pp=fill3(pltmirror(:,1),pltmirror(:,2),pltmirror(:,3),'g');pp.FaceAlpha=0.3;
% 
% pt1=-pltm1*mirror_n1+[0,0,hexpts1(1,3)];
% pt2=pltm1*mirror_n1+[0,0,hexpts1(1,3)];
% pp=line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',[0,1,0,0.3],'LineWidth',7)
% %pp.FaceAlpha=0.3;
% 
% 
% mirror_n1=(hexpts1(2,:)+hexpts1(3,:))/2-(hexpts1(5,:)+hexpts1(6,:))/2;
% mirror_n1=mirror_n1/norm(mirror_n1);
% mirror_n2=[0,0,1];
% 
% pltm1=4.;
% pltm2=1.0;
% pltmirror=zeros(4,3);
% pltmirror(1,:)=-pltm1*mirror_n1-pltm2*mirror_n2;
% pltmirror(2,:)=pltm1*mirror_n1-pltm2*mirror_n2;
% pltmirror(3,:)=pltm1*mirror_n1+pltm2*mirror_n2;
% pltmirror(4,:)=-pltm1*mirror_n1+pltm2*mirror_n2;
% pltmirror(:,3)=pltmirror(:,3)+hexpts1(1,3);
% %pp=fill3(pltmirror(:,1),pltmirror(:,2),pltmirror(:,3),'g');pp.FaceAlpha=0.3;
% 
% pt1=-pltm1*mirror_n1+[0,0,hexpts1(1,3)];
% pt2=pltm1*mirror_n1+[0,0,hexpts1(1,3)];
% pp=line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',[0,1,0,0.3],'LineWidth',7)
% %pp.FaceAlpha=0.3;
% 
% 
% 
% 
% 
% 
% sig_x=[0,1;1,0];
% sig_y=[0,-i;i,0];
% sig_z=[1,0;0,-1];
% sig_0=eye(2);
% 
% gen_spinh=@(nnvec) (sig_x*nnvec(1)+sig_y*nnvec(2)+sig_z*nnvec(3))/norm(nnvec);
% 
% eva_spinS =@(state_wf) [state_wf'*sig_x*state_wf,state_wf'*sig_y*state_wf,state_wf'*sig_z*state_wf];
% 
% 
% 
% newnvec_tetra=Rh_posC(3,:)-tetra1C;
% newnvec_tetra=newnvec_tetra/norm(newnvec_tetra);
% 
% 
% 
% 
% 
% tripts=Rh_posC(1:3,1:3);
% tripts_xaxis=newx_axis(1:3,1:3);
% %quiver3(tripts(:,1),tripts(:,2),tripts(:,3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% for ind=1:3
%     %plt_porbital(tripts(ind,:),tripts_xaxis(ind,:));
% end
% 
% 
% 
% 
% 
% shiftvec=a2;
% tripts=Rh_posC([1,2,4],1:3);
% %tripts_xaxis=newx_axis([1,2,4],1:3);
% 
% newnvec=Rh_posC(3,:)-tetra1_C;
% newnvec=newnvec/norm(newnvec);
% newnvecx=[1,0,0];
% newnvecy=cross(newnvec,newnvecx);
% tripts_xaxis=zeros(3,3);
% tripts_xaxis(1,:)=0.5*newnvecx-0.5*sqrt(3)*newnvecy;
% tripts_xaxis(2,:)=0.5*newnvecx+0.5*sqrt(3)*newnvecy;
% tripts_xaxis(3,:)=[-1,0,0];
% 
% tripts_zaxis=zeros(3,3);
% for ind=1:3
%     tripts_zaxis(ind,:)=tripts(ind,:)-tetra1_C;
% 
% end
% 
% for ind=[1,2]
%     %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
%     plt_D0_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);
% 
% end
% 
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% 
% 
% R3mat=zeros(3);
% R3mat(3,3)=1;
% R3mat(1,1)=-0.5;
% R3mat(2,2)=-0.5;
% R3mat(1,2)=0.5*sqrt(3);
% R3mat(2,1)=-0.5*sqrt(3);
% 
% 
% tripts_xaxis=(R3mat*tripts_xaxis')';
% shiftvec=a1;
% tripts=Rh_posC([3,1,4],1:3);
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% 
% tripts_zaxis=zeros(3,3);
% for ind=1:3
%     tripts_zaxis(ind,:)=tripts(ind,:)-tetra1_C;
% 
% end
% 
% 
% 
% for ind=1:3
%     %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
% end
% 
% for ind=[1,2]
%     %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
%     plt_D0_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);
% 
% end
% 
% 
% tripts_xaxis=(R3mat*tripts_xaxis')';
% shiftvec=a3;
% tripts=Rh_posC([2,3,4],1:3);
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% for ind=1:3
%     %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
% end
% 
% 
% tripts_zaxis=zeros(3,3);
% for ind=1:3
%     tripts_zaxis(ind,:)=tripts(ind,:)-tetra1_C;
% 
% end
% 
% for ind=[1,2]
%     %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
%     plt_D0_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);
% 
% end
% 
% 
% %plt_D2_orbital([0,0,0],[0,0,0]);
% 
% 'add light'
% 
% %set(gca, 'FontSize', 10, 'XColor', 'black', 'YColor', 'black', 'ZColor', 'black', 'LineWidth', 1);
% %axis equal
% %lightangle(-45,30);
% %campos([-1.0926  -1.7493  0.9397]);
% 
% %camlight('headlight')
% 
% 
% 
% 
% 
% 
% 
% % plt_D0_orbital([0,0,1],[1,0,0],[0,1,1],-1);
% 
% %view([1,1,0.5])
% 
% %view([139,22])
% view([0,0,1])
% 
% set(gca,'visible','off')
% 
% %camlight('headlight')
% 
% hold off;
% 
% 
% box on;
% axis equal;
% 
% 
% 
% 
% [az,el] = view



%%


%exportgraphics(gcf,'Pyrochlore_dz2_loopstate.eps','ContentType','vector','BackgroundColor','none')

figure(1)
clf

hold on;


linecc= [.5 .5 .5];
%plot_tetra(Rh_posC,0*a1,linecc)
plot_tetra_v4(Rh_posC,a2,linecc,2)


shiftvec=a2;
tripts=Rh_posC([1,2,4],1:3);
%tripts_xaxis=newx_axis([1,2,4],1:3);

newnvec=Rh_posC(3,:)-tetra1_C;
newnvec=newnvec/norm(newnvec);
newnvecx=[1,0,0];
newnvecy=cross(newnvec,newnvecx);
tripts_xaxis=zeros(3,3);
tripts_xaxis(1,:)=0.5*newnvecx-0.5*sqrt(3)*newnvecy;
tripts_xaxis(2,:)=0.5*newnvecx+0.5*sqrt(3)*newnvecy;
tripts_xaxis(3,:)=[-1,0,0];

tripts_zaxis=zeros(3,3);
for ind=1:3
    tripts_zaxis(ind,:)=tripts(ind,:)-tetra1_C;

end

for ind=[1,2]
    %plt_porbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:));
    plt_D0_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);

end
ind=3
plt_D0_orbital_v3(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),tripts_zaxis(ind,:),(-1)^ind);


pltmirror=[0,-1,-1.3;
    0,-1,0.6;
    0,1,0.6;
    0,1,-1.3]*2.5;

pltmirror(:,1)=pltmirror(:,1)+tetra1_C(1)+shiftvec(1);
pltmirror(:,2)=pltmirror(:,2)+tetra1_C(2)+shiftvec(2);
pltmirror(:,3)=pltmirror(:,3)+tetra1_C(3)+shiftvec(3);


pp=fill3(pltmirror(:,1),pltmirror(:,2),pltmirror(:,3),'g');pp.FaceAlpha=0.3; % 0.3




pt1=tripts(1,:)+shiftvec;
pt2=tripts(2,:)+shiftvec;
line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color','cyan','LineWidth',7)

%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')


R3mat=zeros(3);
R3mat(3,3)=1;
R3mat(1,1)=-0.5;
R3mat(2,2)=-0.5;
R3mat(1,2)=0.5*sqrt(3);
R3mat(2,1)=-0.5*sqrt(3);




view([0,0,1])
view([11,8])

set(gca,'visible','off')

%camlight('headlight')

material dull

hold off;


box on;
axis equal;

[az,el] = view

camlight('headlight')




%%


set(gca,'visible','off')
set(gcf, 'color', 'none');
set(gca, 'color', 'none');


exportgraphics(gcf,'Fig4c-pzCLS-twosite.eps','ContentType','image','BackgroundColor','none','Resolution',400)
% 'ContentType','vector'


% 
% 
% set(gca,'visible','off')
% set(gcf, 'color', 'none');
% set(gca, 'color', 'none');
% 
% 
% exportgraphics(gcf,'Fig4c-twosite-test2.eps','ContentType','vector','BackgroundColor','none')

